Excess mortality: 1877-2020 monthly BfS data & regression approach
Data
Monthly number of deaths, 1877-2020.
Prepared in 02.Rmd.
bfs_web_deaths_monthly <- read_rds("data/BfS/bfs_web_deaths_monthly.Rds") %>%
filter(type == "monthly") %>%
select(-type, -deaths_year) %>%
rename(deaths = deaths_month)Years
year_from <- min(bfs_web_deaths_monthly$Year)
year_reg <- year_from + 10
year_max <- 2020Pandemics
flu_years_hc <- c(1890, 1918, 2020)
flu_years_hc_affected <- c(seq(1890 + 1, 1890 + 10),
seq(1918 + 1, 1918 + 10))
flu_years <- c(1890, 1918, 1957, 1968, 1977, 2009, 2020)
flu_years_affected <- c(flu_years_hc_affected,
seq(1957 + 1, 1957 + 10),
seq(1968 + 1, 1968 + 10),
seq(1977 + 1, 1977 + 10),
seq(2009 + 1, 2009 + 10))Pops
Data till 2019 - taking mean from jan & dec values.
bfs_web_pop_yearly <- read_rds("data/BfS/bfs_web_pop_yearly.Rds") %>%
filter(Year >= min(bfs_web_deaths_monthly$Year)) %>%
mutate(pop = round((pop_jan + pop_dec) / 2)) %>%
select(-pop_jan, -pop_dec)2020 simply predicted using pops from last 10 years
(pop_2020_pred <- predict(
lm(pop ~ Year,
data = bfs_web_pop_yearly %>% filter(Year >= 2010)),
new = tibble(Year = 2020))) 1
8693850
Together
bfs_web_deaths_monthly %<>%
left_join(bfs_web_pop_yearly) %>%
mutate(flu_year_hc = ifelse(Year %in% flu_years_hc, 1, 0),
flu_year_hc_affected = ifelse(Year %in% flu_years_hc_affected, 1, 0),
deaths_flu_year_hc = ifelse(Year %in% flu_years_hc, NA, deaths),
flu_year = ifelse(Year %in% flu_years, 1, 0),
flu_year_affected = ifelse(Year %in% flu_years_affected, 1, 0),
deaths_flu_year = ifelse(Year %in% flu_years, NA, deaths),
pop = ifelse(Year == 2020, pop_2020_pred, pop),
si_one = sin(2*pi*Month/12),
si_two = sin(4*pi*Month/12),
co_one = cos(2*pi*Month/12),
co_two = cos(4*pi*Month/12))
rm(bfs_web_pop_yearly)4 harmonic variables
bfs_web_deaths_monthly %>%
filter(Year >= 2018) %>%
select(si_one, co_one) %>%
mutate(row = row_number()) %>%
pivot_longer(-row) %>%
ggplot() +
geom_line(aes(row, value, color = name)) +
xlab("") +
theme_minimal()bfs_web_deaths_monthly %>%
filter(Year >= 2018) %>%
select(si_two, co_two) %>%
mutate(row = row_number()) %>%
pivot_longer(-row) %>%
ggplot() +
geom_line(aes(row, value, color = name)) +
xlab("") +
theme_minimal()Three pandemics
With denominator
5- & 10-year medians
bfs_web_deaths_monthly %<>%
arrange(Month, Year) %>%
mutate(median = slide_dbl(deaths,
stats::median, na.rm = TRUE,
.before = 10, .after = -1, .complete = TRUE),
median_flu_year_hc = slide_dbl(deaths_flu_year_hc,
stats::median, na.rm = TRUE,
.before = 10, .after = -1, .complete = TRUE),
median_flu_year = slide_dbl(deaths_flu_year,
stats::median, na.rm = TRUE,
.before = 10, .after = -1, .complete = TRUE),
) %>%
arrange(Year, Month)Impact of 1957 on 1958
Seasonality
dcmp <- bfs_web_deaths_monthly %>%
mutate(Time = row_number()) %>%
tsibble::as_tsibble(index = Time) %>%
fabletools::model(feasts::STL((deaths/pop)*100000 ~ season(window = Inf)))
# fabletools::components(dcmp)
fabletools::components(dcmp) %>%
autoplot()rm(dcmp)Outcome
Marcel’s regression approach
CI via bootstrap
Method from https://statcompute.wordpress.com/2015/12/20/prediction-intervals-for-poisson-regression/
p_load(doParallel, foreach)
registerDoParallel(cores = 4)
boot_pi <- function(model, pdata, n, p) {
odata <- model$data
lp <- (1 - p) / 2
up <- 1 - lp
seeds <- round(runif(n, 1, 1000), 0)
boot_y <- foreach(i = 1:n, .combine = rbind) %dopar% {
set.seed(seeds[i])
bdata <- odata[sample(seq(nrow(odata)), size = nrow(odata), replace = TRUE), ]
bpred <- predict(update(model, data = bdata), type = "response", newdata = pdata)
rpois(length(bpred), lambda = bpred)
}
boot_ci <- t(apply(boot_y, 2, quantile, c(lp, up)))
return(data.frame(pred = predict(model, newdata = pdata, type = "response"),
lower = boot_ci[, 1], upper = boot_ci[, 2]))
}Looping over years
Analyses using 10 previous years.
With three different inclusions depending on pandemic / flu years
for (YEAR in year_reg:year_max){
print(paste("Analysing year", YEAR))
reg_data <- bfs_web_deaths_monthly %>%
filter(Year >= YEAR - 10 & Year < YEAR)
pred_data <- bfs_web_deaths_monthly %>%
filter(Year == YEAR) %>%
mutate(
fit = NA_real_,
lpi = NA_real_,
upi = NA_real_,
fit_flu_hc = NA_real_,
lpi_flu_hc = NA_real_,
upi_flu_hc = NA_real_,
fit_flu = NA_real_,
lpi_flu = NA_real_,
upi_flu = NA_real_)
# summary(reg_data$Year)
# summary(pred_data$Year)
# all data
summary(monthly <- glm(deaths ~ Year +
si_one + si_two + co_one + co_two,
offset = log(pop),
data = reg_data,
family = "poisson"))
predict <- boot_pi(monthly, pred_data, 1000, 0.99)
pred_data %<>%
mutate(
fit = predict$pred,
lpi = predict$lower,
upi = predict$upper,
fit_flu_hc = predict$pred,
lpi_flu_hc = predict$lower,
upi_flu_hc = predict$upper,
fit_flu = predict$pred,
lpi_flu = predict$lower,
upi_flu = predict$upper
)
# flu hc excluded
if (YEAR %in% flu_years_hc_affected) {
print(paste(" >> Extra analyses for hc flu year exclusions"))
summary(monthly_flu_hc <- glm(deaths_flu_year_hc ~ Year +
si_one + si_two + co_one + co_two,
offset = log(pop),
data = reg_data,
family = "poisson"))
predict_flu_hc <- boot_pi(monthly_flu_hc, pred_data, 1000, 0.99)
pred_data %<>%
mutate(
fit_flu_hc = predict_flu_hc$pred,
lpi_flu_hc = predict_flu_hc$lower,
upi_flu_hc = predict_flu_hc$upper
)
}
# all flu excluded
if (YEAR %in% flu_years_affected) {
print(paste(" >> Extra analyses for flu year exclusions"))
summary(monthly_flu <- glm(deaths_flu_year ~ Year +
si_one + si_two + co_one + co_two,
offset = log(pop),
data = reg_data,
family = "poisson"))
predict_flu <- boot_pi(monthly_flu, pred_data, 1000, 0.99)
pred_data %<>%
mutate(
fit_flu = predict_flu$pred,
lpi_flu = predict_flu$lower,
upi_flu = predict_flu$upper
)
}
if (YEAR == year_reg) {
expected_deaths_monthly <- pred_data
} else {
expected_deaths_monthly <- bind_rows(expected_deaths_monthly, pred_data)
}
}[1] "Analysing year 1887"
[1] "Analysing year 1888"
[1] "Analysing year 1889"
[1] "Analysing year 1890"
[1] "Analysing year 1891"
[1] " >> Extra analyses for hc flu year exclusions"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1892"
[1] " >> Extra analyses for hc flu year exclusions"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1893"
[1] " >> Extra analyses for hc flu year exclusions"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1894"
[1] " >> Extra analyses for hc flu year exclusions"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1895"
[1] " >> Extra analyses for hc flu year exclusions"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1896"
[1] " >> Extra analyses for hc flu year exclusions"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1897"
[1] " >> Extra analyses for hc flu year exclusions"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1898"
[1] " >> Extra analyses for hc flu year exclusions"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1899"
[1] " >> Extra analyses for hc flu year exclusions"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1900"
[1] " >> Extra analyses for hc flu year exclusions"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1901"
[1] "Analysing year 1902"
[1] "Analysing year 1903"
[1] "Analysing year 1904"
[1] "Analysing year 1905"
[1] "Analysing year 1906"
[1] "Analysing year 1907"
[1] "Analysing year 1908"
[1] "Analysing year 1909"
[1] "Analysing year 1910"
[1] "Analysing year 1911"
[1] "Analysing year 1912"
[1] "Analysing year 1913"
[1] "Analysing year 1914"
[1] "Analysing year 1915"
[1] "Analysing year 1916"
[1] "Analysing year 1917"
[1] "Analysing year 1918"
[1] "Analysing year 1919"
[1] " >> Extra analyses for hc flu year exclusions"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1920"
[1] " >> Extra analyses for hc flu year exclusions"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1921"
[1] " >> Extra analyses for hc flu year exclusions"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1922"
[1] " >> Extra analyses for hc flu year exclusions"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1923"
[1] " >> Extra analyses for hc flu year exclusions"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1924"
[1] " >> Extra analyses for hc flu year exclusions"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1925"
[1] " >> Extra analyses for hc flu year exclusions"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1926"
[1] " >> Extra analyses for hc flu year exclusions"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1927"
[1] " >> Extra analyses for hc flu year exclusions"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1928"
[1] " >> Extra analyses for hc flu year exclusions"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1929"
[1] "Analysing year 1930"
[1] "Analysing year 1931"
[1] "Analysing year 1932"
[1] "Analysing year 1933"
[1] "Analysing year 1934"
[1] "Analysing year 1935"
[1] "Analysing year 1936"
[1] "Analysing year 1937"
[1] "Analysing year 1938"
[1] "Analysing year 1939"
[1] "Analysing year 1940"
[1] "Analysing year 1941"
[1] "Analysing year 1942"
[1] "Analysing year 1943"
[1] "Analysing year 1944"
[1] "Analysing year 1945"
[1] "Analysing year 1946"
[1] "Analysing year 1947"
[1] "Analysing year 1948"
[1] "Analysing year 1949"
[1] "Analysing year 1950"
[1] "Analysing year 1951"
[1] "Analysing year 1952"
[1] "Analysing year 1953"
[1] "Analysing year 1954"
[1] "Analysing year 1955"
[1] "Analysing year 1956"
[1] "Analysing year 1957"
[1] "Analysing year 1958"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1959"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1960"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1961"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1962"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1963"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1964"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1965"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1966"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1967"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1968"
[1] "Analysing year 1969"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1970"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1971"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1972"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1973"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1974"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1975"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1976"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1977"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1978"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1979"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1980"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1981"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1982"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1983"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1984"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1985"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1986"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1987"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 1988"
[1] "Analysing year 1989"
[1] "Analysing year 1990"
[1] "Analysing year 1991"
[1] "Analysing year 1992"
[1] "Analysing year 1993"
[1] "Analysing year 1994"
[1] "Analysing year 1995"
[1] "Analysing year 1996"
[1] "Analysing year 1997"
[1] "Analysing year 1998"
[1] "Analysing year 1999"
[1] "Analysing year 2000"
[1] "Analysing year 2001"
[1] "Analysing year 2002"
[1] "Analysing year 2003"
[1] "Analysing year 2004"
[1] "Analysing year 2005"
[1] "Analysing year 2006"
[1] "Analysing year 2007"
[1] "Analysing year 2008"
[1] "Analysing year 2009"
[1] "Analysing year 2010"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 2011"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 2012"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 2013"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 2014"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 2015"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 2016"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 2017"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 2018"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 2019"
[1] " >> Extra analyses for flu year exclusions"
[1] "Analysing year 2020"
rm(predict, predict_flu_hc, predict_flu, pred_data, reg_data, YEAR)
write_rds(expected_deaths_monthly, "data/expected_deaths_monthly.Rds")Two predictions
1890
With extremes
Without extremes
1918
With extremes
Without extremes
1957
With extremes
Without extremes
2020
With extremes
Without extremes
Alt Viz
Animate
Stacked
Diffs
Computing Environment
R version 4.0.4 (2021-02-15) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows 10 x64 (build 18363) Matrix products: default attached base packages: [1] parallel stats graphics grDevices utils datasets methods [8] base other attached packages: [1] iterators_1.0.13 slider_0.2.1 magrittr_2.0.1 scales_1.1.1 [5] lubridate_1.7.10 forcats_0.5.1 stringr_1.4.0 dplyr_1.0.5 [9] purrr_0.3.4 readr_1.4.0 tidyr_1.1.3 tibble_3.1.1 [13] ggplot2_3.3.3 tidyverse_1.3.1 pacman_0.5.1To cite R in publication use:
R Core Team (2021). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/.